require(stats)
require(graphics)
#library("ggplot2")
require(ggplot2)
#library("grDevices")
require(grDevices)
#library("gridExtra")
require(gridExtra)
#library("grid")
require(grid)
#install.packages("Hmisc")
require(Hmisc)

# Please make sure to change the path to where the appropriate dataset is located on your computer

# FULL RESTRICT
fr.k1.data <- read.csv("PNAS_VARILL_FR_K1.csv")
# STUDENT
st.k1.data <- read.csv("PNAS_VARILL_ST_K1.csv")
# FAMILY
fam.k1.data <- read.csv("PNAS_VARILL_FAM_K1.csv")
#LOW SKILLED
ls.k1.data <- read.csv("PNAS_VARILL_LS_K1.csv")
#HIGH SKILLED
hs.k1.data <- read.csv("PNAS_VARILL_HS_K1.csv")
#BASELINE
open.k1.data <- read.csv("PNAS_VARILL_O_K1.csv")
#FREE MOVEMENT
no.k1.data <- read.csv("PNAS_VARILL_NO_K1.csv")

sponsor.require <- 1
high.skilled.p <-  2
low.skilled.p <-  3
p.tourist.visa <- 4
p.fakedocs <-  5
no.econ.stud <- 6
p.all.illegal <- 7
k <- 8
my.seed <- 9

# Index columns

legaltourist.attempts <- c(10:29)

ability.legal.cols <- c(30:90)
ability.fakedocs.cols <- c(91:151)
ability.work.cols <- c(152:212)

legal.cols <- c(274:334) # legal
fi.cols <- c(335:395) # illegal
sl.cols <- c(396:456) # semi-legal
stud.mig <- c(457:517) # student
hs.mig <- c(518:578) # high-skilled
ls.mig <- c(579:639) #low-skilled
fam.mig <- c(640:700) #family

# FROM INPUT DATA IN SEPARATE FILE
num.agents <-  1166 # all agents
aspire.enough <- 507 # agents who aspire to migrate

table(fr.k1.data[,p.all.illegal])
table(st.k1.data[,p.all.illegal])
table(fam.k1.data[,p.all.illegal])
table(ls.k1.data[,p.all.illegal])
table(hs.k1.data[,p.all.illegal])
table(open.k1.data[,p.all.illegal])
table(no.k1.data[,p.all.illegal])

select.policy.level <-  0.3 #1-0.7

# LIST OF ALL DATASETS
df.list <- list(no.k1.data,fr.k1.data,fam.k1.data,ls.k1.data,open.k1.data,hs.k1.data,st.k1.data)

# USE ENFORCEMENT SETTING FOR FIGURE 2
df.list <- lapply(df.list, function(x) {
  x <- x[which(x[,p.all.illegal] == select.policy.level),]
  x} )

# MAKE ALL DATASETS BE OF EQUAL LENGTH FOR SUBTRACTION
df.list <- lapply(df.list, function(x) {
  x <- x[1:982,]
  x} )

names.df.list <-  Cs(no.k1.data,fr.k1.data,fam.k1.data,ls.k1.data,open.k1.data,hs.k1.data,st.k1.data)

names(df.list) <- names.df.list

# YEAR AT WHICH WE COMPUTE
time.point <- 20

for(i in 1:length(df.list))
{
  df.list[[i]]$legal.raw <- (df.list[[i]][,legal.cols[time.point]])
  df.list[[i]]$legal <-  df.list[[i]]$legal.raw / aspire.enough
  df.list[[i]]$fill <- (df.list[[i]][,fi.cols[time.point]]) #illegal
  df.list[[i]]$semi <- (df.list[[i]][,sl.cols[time.point]]) #semi-legal
  df.list[[i]]$all.illegal.raw <- (df.list[[i]]$fill + df.list[[i]]$semi)
  df.list[[i]]$all.illegal <- df.list[[i]]$all.illegal.raw / aspire.enough
  df.list[[i]]$nm <- 1 - (df.list[[i]]$legal + df.list[[i]]$all.illegal) #non-migrants
  df.list[[i]]$all.migs.raw <- df.list[[i]]$legal.raw + df.list[[i]]$all.illegal.raw
}

# Subtract models from baseline (open)

b.i <- grep("open", names.df.list)

for (i in 1:7)
{
  df.list[[i]]$legal.adj <- df.list[[i]]$legal - df.list[[b.i[1]]]$legal
  df.list[[i]]$nm.adj  <- df.list[[i]]$nm - df.list[[b.i[1]]]$nm
  df.list[[i]]$illegal.adj  <- df.list[[i]]$all.illegal - df.list[[b.i[1]]]$all.illegal
  df.list[[i]]$legal.migs.adj <- (df.list[[i]]$legal.raw - df.list[[b.i[1]]]$legal.raw) / df.list[[i]]$all.migs.raw
  df.list[[i]]$illegal.migs.adj  <- (df.list[[i]]$all.illegal.raw - df.list[[b.i[1]]]$all.illegal.raw) / df.list[[i]]$all.migs.raw
}

k1 <- matrix(nrow = 7, ncol = 3)

k1.mig <- matrix(nrow = 7, ncol = 2)

for (i in 1:7)
{
  k1[i,1] <- mean(df.list[[i]]$illegal.adj) 
  k1[i,2] <- mean(df.list[[i]]$nm.adj) 
  k1[i,3] <- mean(df.list[[i]]$legal.adj) 
  k1.mig[i,1] <- mean(df.list[[i]]$illegal.migs.adj) 
  k1.mig[i,2] <- mean(df.list[[i]]$legal.migs.adj)
}

k1 <- as.table(k1[-b.i[1],])
k1 <-  k1 * 100 # percentages

k1.mig <- as.table(k1.mig[-b.i[1],])
k1.mig <- k1.mig * 100 #percentages

tnames <- c("Free\nMovement","Closed","Family","Low-Skilled", "High-Skilled","Student")
tnames.empty <- c("","","","","","")


pdf("Fig2.pdf",width=6,height=6)
par(pty="s", mar=c(3, 7, 2, 2), lwd = 0.2)
barplot(t(k1[,c(1,2)]), horiz=T, las=1, col=c("firebrick1","black"), names.arg = tnames,
        border = T,cex.names=1.2, cex.axis = 1.2, font.axis = 1.2,
        cex.lab = 1.2, xlab = "% of Aspiring Migrants", xlim = c(-60,60), axes = T, 
        legend.text=c("Unauthorized","Non Migrants","Legal"), 
        args.legend= list(fill=c("firebrick1","black","lightblue")))

barplot(t(k1[,3]), horiz=T, las=1, col=c("lightblue"), names.arg = tnames.empty,
        border = T,cex.names=1.2, cex.axis = 1.2, font.axis = 1.2, axes = F, add = T)

dev.off()
